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Granular matter is comprised of a large number of particles whose collective behavior determines 
macroscopic properties such as flow and mechanical strength [l], Q. A comprehensive theory of the 
properties of granular matter, therefore, requires a statistical framework. In molecular matter, 
equilibrium statistical mechanics [3], which is founded on the principle of conservation of energy, 
provides this framework. Grains, however, are small but macroscopic objects whose interactions are 
dissipative since energy can be lost through excitations of the internal degrees of freedom. In this 
work, we construct a statistical framework for static, mechanically stable packings of grains, which 
parallels that of equilibrium statistical mechanics but with conservation of energy replaced by the 
conservation of a function related to the mechanical stress tensor. Our analysis demonstrates the 
existence of a state function that has all the attributes of entropy [3J. In particular, maximizing this 
state function leads to a well-defined granular temperature for these systems. Predictions of the 
ensemble are verified against simulated packings of frictionless, deformable disks. Our demonstration 
that a statistical ensemble can be constructed through the identification of conserved quantities other 
than energy is a new approach that is expected to open up avenues for statistical descriptions of 
other non-equilibrium systems. 



Phases and phase transitions are manifestations of the 
collective behavior of many-particle systems. Phases are 
clearly delineated in atomic and molecular systems where 
solids resist shear but liquids do not. In granular mat- 
ter, however, phases lose some of their distinction [![: a 
pile of sand or coffee beans behaves as a solid but can 
break down in avalanches with only the surface layer 
flowing. Ordinary matter is a statistical ensemble of 
molecules with energy conserving interactions. In such 
systems, equilibrium statistical mechanics connects the 
microscopic properties to the observable, macroscopic be- 
havior through probability distributions characterized by 
parameters such as temperature, and density 0]. Gran- 
ular materials differ from their molecular counterparts 
in two important respects; (a) the grains are themselves 
macroscopic and lose energy every time they interact and 
(b) because of their size, gravitational effects swamp ther- 
mal fluctuations and changes are effected through me- 
chanical acts such as shaking and tumbling. The dissipa- 
tive nature of granular materials makes them explicitly 
non-equilibrium: energy has to be supplied to the sys- 
tem to maintain a steady state. If the energy supply 
is turned off, granular materials relax to a static, me- 
chanically stable packing (a blocked state) through the 
dissipation of energy Response and phase behavior of 
weakly-driven granular systems are determined by the 
excitations of these blocked states. Prediction of these 
collective phenomena entails making a statistical connec- 
tion between the grain-level and macroscopic properties 
of blocked states since, given a fixed set of macroscopic 
parameters, e.g. volume and number of grains, many 
such states are possible . 

In this work, we construct a statistical framework 
for predicting probability distributions of blocked states, 
which applies to frictional and frictionless grains but is 



limited to isotropic ones. The detailed analysis is re- 
stricted to planar packings where the mechanical stress 
tensor is known to have a simple structure 0, 0). The 
fundamental postulates, however, apply to higher dimen- 
sions. 

The paper is organized as follows: (a) through the in- 
troduction of a height field for planar packings, we show 
that the total stress tensor is a conserved quantity; (b) 
we formulate a statistical theory of blocked states which 
parallels that of thermal systems but with the internal 
virialQ, T, playing the role of energy; (c) we postulate 
entropy maximization to introduce a "granular temper- 
ature" and construct a probability distribution involv- 
ing the inverse granular temperature, a ; (d) we check 
this postulate against packings generated through a dy- 
namical algorithm by showing that all parts of a given 
packing have the same a; (e) finally, we use the deduced 
form of a(T) to demonstrate that the entropy is indeed 
the logarithm of the density of states. This last result 
provides strong support for the microcanonical ansatz 
that all states with the same T are equally likely and 
connects to the Edwards hypothesis for incompressible 
grains that all blocked states with the same volume are 
equally likely [7|. 

Every grain in a mechanically stable packing of fric- 
tional or frictionless grains has to satisfy the constraints 
of force and torque balance. For planar packings of grains 
which interact only through contact forces, these con- 
straints can be incorporated through a mapping to a set 
of auxiliary (height) variables defined on the dual net- 
work of the voids surrounded by grains [J Q. Choos- 
ing the center of an arbitrary void as the zero of the 
height field, the height vectors are constructed iteratively 
through: h v = Fij + where ij is the contact traversed 
in going from the center of the void fi to that of v (c/Fig 
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FIG. 1: Height map: A height map with its origin at O in the 
interior of a packing obtained from simulations. The height vectors 
(light arrows) are shown along with the boundary vectors f^i, 
fii2 (dark arrows) entering eq. J7J. The dashed lines denote in- 
tergrain contacts and the ones included in a a for a 8-grain cluster 
occupying the shaded region with area A, are in bold. 
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FIG. 2: Numerical test of the equality of a: Row 1: Results 
from packings with N = 4096, m = 8 and 4> = 0.838-0.844. Row 2: 
Same for cf> = 0.95 and 1.0. First column shows P m ,p(r m ); second 
column is r p , q vs. T m —T' m for an arbitrarily chosen pair (p, q). The 
error bars derive from the spread in r Pi9 obtained for different T m 
and rj„, contributing to a given T m — T^; third column illustrates 
the scaling of P m , p (r m ) according to Eq. [2] 



[T]). Since the i'y's around a grain sum to zero, the map- 
ping of forces to heights is one to one up to an arbitrary 
choice of the origin (O in Fig [T]) . This mapping satis- 
fies Newton's third law and enforces the force-balance 
constraint. The constraint of torque balance imposes a 
divergenceless condition on the height field [J, Q. We use 
this mapping to demonstrate a conservation law related 
to the mechanical stress. 

The height field, which enforces that the stress tensor 
is divergence-free [4], is an exact analog of the vector po- 
tential in clcctromagnetism, which ensures that the mag- 
netic field is divergence- free |9j. As shown in the Methods 
section, the height map leads to an explicit implementa- 
tion of Stokes' theorem Q which proves that the stress 
tensor a a can be expressed as a sum involving only the 
boundary of the area A occupied by the grains. This re- 
sult allows us to define an extensive quantity £ = Ad a 
which is left invariant by local rearrangements. Differ- 
ent packings with the same value of £ can be related 
to each other through processes that do not involve the 
boundary. An example is provided by the "wheel moves" 
in a triangular lattice [Io|. which transposes one force- 
balanced packing into another leaving £ unchanged. 

In two dimensions, the tensor £ has two scalar invari- 
ants and these can be taken to be the trace T, which is 
the internal virialjij], and the determinant, k. For d > 3, 
the loops around grains cannot be defined unambiguously 
and the height map cannot be constructed. The total 
stress tensor, S, is nevertheless conserved under internal 
rearrangements since it can be written as the curl of a 
2 nd rank tensor ll|, and Stokes' theorem can be used to 
represent it as a boundary integral. 

The above analysis demonstrates that for blocked 



states of frictional and frictionless packings of deformable 
disks, T = J2ij dijF.ij (cf Eq. |5J) is a pure boundary 
term. Since packings with different values of T cannot 
be transformed into one another through internal rear- 
rangements not involving the boundary, T becomes the 
analog of energy in thermal systems. It is then rigor- 
ously possible to calculate the density of states 0(r, V, N) 
of all blocked configurations characterized by a given T, 
volume V and number of grains N. Analogous calcu- 
lations, based on the conservation of physical variables 
other than energy, have been performed in certain lattice 
models with rigid constraints An entropy function 
can be defined as S = lnfi(r, V, N) and assumption of 
an entropy maximization principle for granular equilib- 
rium implies equalization of the granular temperature [3j: 
a = dS/dT\v.N- A constant-a canonical ensembleQ fol- 
lows with the probability of finding a blocked state v 
being given by: P v — exp(—aT u )/Z(a,V 1 N). In the 
current work, we ignore the invariant determinant n. It 
can easily be included as S(T, k, V, N) and gives rise to an 
additional factor e~ MK in the probability. For frictionless 
systems, k is directly related to T. 

To test the predictions of our statistical framework, we 
generate sets of mechanically stable packings of a fixed 
number N of deformable circular disks and group them 
according to their packing fraction <p = Nv/V, where 
v is the microscopic grain volume and V is the volume 
of the packing (see Methods section). Each packing is 
characterized by Tat, the total value of T for the iV-grain 
packing and by (z), the average value of the number of 
contacts in the packing. Units of length and interaction 
potential are chosen to make T dimensionless (see Meth- 
ods section). 

The entropy maximization hypothesis implies [3J that 
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all subregions of a grain packing (p), generated through 
some dynamics, have the same granular temperature a p . 
The values of T m for m-grain clusters inside packing p 
should then be distributed according to: P m ,p(Y m ) = 



"(5(r zy — T TO ). By measuring P TOjP at a pair of 



values, r m and T' m , and comparing these to P m ,q for the 
same pair we construct the ratio: 



= loe 



Pm,p (J- *m)Pm,q (F 'rn 
Pm ,p(X ' m ) Pm ,qfi 



= -(a p -a q )(T m -r' m ) 
(1) 

The last equality follows if the packings p and q have 
the same or similar values of <p (see Methods section). 
Packings with N = 4096 in two narrow packing fraction 
ranges are shown in the first column of Figure [2] We find 
that for pairs (p, q) of packings in the low (0.838 — 0.844) 
packing fraction range and in the range (0.95, 1.00), r p<q 
is indeed a linear function of T m — T' m , establishing that 
a is equal for all parts of a grain packing. An example 
of the linear correlation is shown in column 2 of Figure 
[2j The average values of extracted from 

these plots are « 2 x 10 4 and 20, for the low packing 
fractions and the high packing fractions, respectively. A 
further test of the equality of a inside a packing is to 
verify that for a pair (p, q) with negligible differences in 



Pm,,q(J^ni) 6Xp( 5QZ p qY m ) Prn,p(J^m) 



(2) 



if both distributions are normalized to unity. The third 
column in Figure [2] demonstrates the scaling of packings 
with an arbitrary packing chosen as the reference p. 

For packings with tfi ~ 0.84, it is possible to find over- 
lapping distributions, P m ,p{Y m ) over a range of T^r ex- 
tending from ~ 10~ 5 — 10~ 3 . In this regime, we extract 
8ap. q systematically from Eq. [2] by rescaling the same 
packing q with two different reference packings, p, and 
determine a as a function of Tat, up to an arbitrary con- 
stant C. For m > 2, we find that over two orders of 
magnitude in Tat, the results are independent of m and 
described by the functional form: 



a = C + 2N/T 



N 



(3) 



Eq. [3] implies that the granular temperature 1 /a tends 
to for Tn — > 0. The value of C is determined by the 
behavior of a for large Tn . Since we expect I /a — > oo 
as Tn — > oo, we choose C = 0. 

The above numerical tests demonstrate that for all the 
packings studied, over the full range of <fi, different sub- 
regions of a packing have the same value of a. Moreover, 
for r^r < 10~ 3 , a is independent of m and a unique 
function of T^. These results establish a as the granu- 
lar equivalent of inverse temperature. The parallel with 
equilibrium statistical mechanics is further illustrated be- 
low through the verification of the analog of the Boltz- 
man entropy formula^ S cx In Q providing support for 
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FIG. 3: Test of Eq. g) Distributions for N = 4096 and m = 24 
plotted as a function of x = NT m /T N . For <f> = 0.838 to <f> = 0.844, 
distributions are compared to Eq. [4] with a = 2 (solid line). For 
4> = 0.95, with (z) in the range 5.207 ± 0.06, Eq. |4]with a = 5.75 
(solid line) best describes the data. 



the microcanonical ansatz of equal probability of blocked 
states with the same T. 

Since a is defined by the derivative of the entropy with 
respect to T, Eq. [3]implies a specific form of the density of 
states: Q(T, V, N) = e s - T Na , with a = 2, independent 
of <j), for packings with cf> ~ 0.84. Combining this relation 
with the functional form of a(T), it is straightforward to 
show that the distribution of T m in a packing (c/Eq. [8|) 
is a function only of x — NY m /Y n'- 



P m , p (T m ) = P m (x) = Cx ma exp(-ax) 



(4) 



The scaling of the distributions, obtained from multiple 
packings with (f> ~ 0.84, is illustrated in Figure [3] for 
to = 24 and verifies that S is indeed proportional to In £1 
in this regime. 

For packings with <f> > 0.85, we found that the distri- 
butions obtained from different packings no longer scale 
with x — NT m /T]sr but depend explicitly on </>; becoming 
narrower with increasing <fi. We know from the simula- 
tion data and previous work 13] that the average num- 
ber of contacts, (z), is a strongly increasing function of 
the packing fraction (f> and that, in the low packing frac- 
tion regime, (z) — > zi so = 4 (for frictionless disks the 
isostatic packing has 2d contacts where d is the spa- 
tial dimension 3, IH ) . It is also clear from Eq. [4] that 
1 j a is proportional to the width of the rescaled distribu- 
tions. Based on these observations we hypothesize that 
the distributions are still described by Eq. [4] but with 
a = 2 + f((z) — z iso ) and /(0) = 0. A test of this hypoth- 
esis is presented in Figured! For packings with tfi = 0.95, 
the measured (z) « 5.2 and the best fit to Eq. |4] yields 
a = 5.75. For smaller systems, with TV = 1024, we per- 
formed numerical fits to the distributions over a wide 
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FIG. 4: Numerical test of the model for a: Fitting parameter 
a obtained from distributions of N = 1024 grain-packings with 
different values of the average number of contacts, (z). 



range of (z). Figure [4] shows the numerical values of a 
for a discrete set of (z) values and for 8 < m < 128. It is 
clear that a increases with (z) and approaches a constant 
for (z) — > 4. 

The results presented above show that the granular 
temperature 1/a of a A^-grain packing is determined by 
(z) and Tpf: a ~ Na((z))/T^. This relationship is the 
analog of an equation of state relating energy, tempera- 
ture and density in a thermal system. Since a approaches 
a constant as Tjv — ^ and (z) — > Zi SO , in this limit, the 
density of states is, to a good approximation, only a func- 
tion of T: 



n(v,v,N) ~ r 2iV oc f 



7 2N 



(5) 



generic dynamics with no a priori information about the 
ensemble and, therefore we expect our results to be ap- 
plicable in general. 

The statistical framework presented in this paper leads 
to a natural phase space for jamming[18] with the axes 
being thermodynamic temperature (T) , packing fraction 
(</>), and granular temperature 1/a. Simulations [T^ and 
experiments[19] have indicated the existence of a T = 
critical point (Point J) on the <p axis. The a ensemble 
provides a framework for exploring phases, phase transi- 
tions and critical points on the T — plane. A prelimi- 
nary attempt, based on the present framework but with- 
out any rigorous justification of its fundamental premise, 
has been made by two of us[8[ and has led to a field the- 
ory of Point J. One of the crucial aspects of this field 
theory is the coupling between (z) and T. Our current 
analysis shows that this coupling is a consequence of the 
nature of the density of states and that the geometry and 
stress decouple only in the infintely-rigid-grain or the in- 
compressible limit (c/Eq. [5]); i.e. at Point J. 

The full richness of granular statistical mechanics, 
which is founded on a conserved tensor not a scalar, 
can be explored only in systems with friction. Our tests 
of the statistical ensemble so far have been confined to 
numerical simulations of frictionless disks because the 
verification required extensive statistical sampling. We 
are in the process of extending the verification of the 
statistical ensemble to experimental packings of grains 
with friction[2Cj|. 

Methods 



Height Map The microscopic stress tensor for grain i is defined 
as 211 



Here F 



. ,.j is the sum of the magnitude of all 
the contact forces and we have used the fact that as 



(6) 



r — > 0, the compression of the grains approach zero and, where the dy are the vectors connecting the grain center to its 



in this limit, T = J2ijdijFij oc F. The statistical en- 
semble of packings in this incompressible limit should 
be identical to that of infinitely rigid grains since, as 
we have demonstrated, the ensemble is determined by 
r^v- The conservation of T reduces to a conservation of 
F for infinitely rigid grains and we recover the Edwards- 
microcanonical ensemble based on total volume and total 
external force @, EH, 17 1. Such an ensemble would lead to 
a density of states of the form ([5]) . 

In summary, we have shown that it is possible to 
construct a statistical ensemble for mechanically stable 
(blocked) states of frictional and frictionless granular 
packings. Numerical tests demonstrate the existence of a 
granular temperature which is given by the logarithmic 
derivative of the density of states with respect to T. The 
tests also establish that the density of states is history- 
independent since different paths leading to the same 
generate the same distributions. Our analysis was per- 
formed on data obtained from a specific though fairly 



contact points, the Fij are the contact forces and Zi is the number 
of contacts for grain i. The height map can be used to rewrite 
&t as[j: &i = E^Lifti +f )1 2)h^, where the r M i and r M 2 are the 
vectors linking the two contact points associated with the the void fj, 
(c/Fig[T]l. The grain area is defined as the area enclosed by the re- 
verters, and these tesselate the plane (c/FigflJ. The macroscopic 
stress tensor a is related to its single grain counterpart through 
, where the summation is over a connected cluster 
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of grains and the area A is the combined grain area of all grains in 
the cluster. In the summation over grains in a a, the terms related 
to voids fully inside the area A add up to zero and, as illustrated 
in Fig[T] we are left with a sum over the voids at the boundary of 
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Simulations We generate mechanically stable packings of bidis- 
perse deformablc disks that interact via purely repulsive linear 
spring interactions for systems with N = 1024 and N = 4096 
disks. The mixtures are 50-50 by number and the diameter 
ratio between large and small disks is 1.4. For each of these 
sizes we study several packing fractions from around random 
close packing <j> = 0.84 to more than 20% above this value. To 
create the packings, we initialize the disks with random initial 
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conditions at a specified packing fraction and then implement 
conjugate gradient energy minimization to find the nearest local 
energy minimum 1 13| . Near random close packing, the packings 
typically have some grains with no contacts ("rattlers"). In 
all of our analysis, these rattlers have been excluded. In the 
simulations, lengths are measured in units of the large-particle 
diameter and ener gy is measured in units of the characteristic 
interaction strength 1 13H . This renders T dimensionless and the con- 
tact force is numerically equal to the magnitude of the disk overlap. 

Testing equality of granular temperature We divide 
a packing p into subregions containing m grains and con- 
struct the probability distribution P m (r m ) by histogramming 
r m = 5ZijCm r ijFij is summed over all contacts in the subregion 
(including boundary contacts). 

The prediction of the canonical a-ensemble is that this distribution 
should have the form 

p m , p (r m ) = ^e- Q * r "5(iv-r m ) 

V 

= Q(r rn ,Vm,m)exp(-a p r rn )/Zm(ap) (8) 

where Q(T m ,v m ,m) is the density of states corresponding to 
v m , the volume of the subregion and Z m (a p ) is the partition 
function^. In practice, we evaluate P m ,p(r m ) by combining all 
m-grain regions irrespective of the value of v m . The appropriate 
distribution to compare the numerical data to is, therefore, 

P m , p (T m ) = Uff/ (r m ,V p ,N)e x p(-a p r m )/Z m (a p ) (9) 

where Q e j f depends on the volume V p of the packing p but not on 
the local geometry. For a given pair of packings with with the same 
N and packing fraction differences small enough such that the V 
dependence of O e ** can be ignored, 

r V,q = ( p Tri~YP TfT^ J - "( a P ~ a lK V rn ~ T m ) (10) 

since the factors of fl e ff from the numerator and denominator 
cancel out. 
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